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ABSTRACT 
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The tight correlation between galaxy bulges and their central black hole 
masses likely emerges in a phase of rapid collapse and starburst at high redshift, 
due to the balance of gravity on gas with the feedback force from starbursts and 
the wind from the black hole; the average gravity on per unit mass of gas is 
' ~ 2 x lCU 10 m sec -2 during the star burst phase. This level of gravity could come 

from the real r _1 cusps of Cold Dark Matter (CDM) halos, but the predicted 
^ ■ gravity would have a large scatter due to dependence on cosmological parame- 

ters and formation histories. Better agreement is found with the gravity from the 
scalar field in some co-variant versions of MOND, which can create the mirage 
of a Newtonian effective dark halo of density Ur~ l near the center, where the 
characteristic surface density II = 13Oa _1 M pc -2 and a is a fundamental con- 
C<") ■ stant of order unity fixed by the Lagrangian of the co-variant theory if neglecting 

^ ■ environmental effects. 

We show with a toy analytical model and a hydrodynamical simulation that a 
constant background gravity due to MOND/TeVeS scalar field implies a critical 
CSj ■ pressure synchronizing starbursts and the formation of galaxy bulges and ellipti- 

cals. A universal threshold for the formation of the brightest regions of galaxies in 
a MONDian universe suggests that the central BHs, bulges and ellipticals would 
respect tight correlations like the Mb u i ge — Mbh — a relations. In general MOND 
^ tends to produce tight correlations in galaxy properties because its effective halo 

has less freedom and scatter than CDM halos. 
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1. Correlated formation of Black Hole and Bulges in CDM and MOND 

While appearing in wide range of shapes, sizes and luminosities, galaxies have very 
regular properties. E.g., the terminal rotation speed V C i r of a spiral galaxy is tightly cor- 
related with its total baryonic mass M, following a simple TFMM power-law Vj r /M ~ 
O.O2(kms _1 ) _4 M _1 , a formula proposed by Tully & Fisher (1977) for high-surface bright- 
ness galaxies, and generalized by Milgrom (1983) and tested by McGaugh (2005) for gas-rich 
low-surface brightness galaxies. There is no evidence for any significant scatter, and it seems 
to apply independent of galaxy formation history (Gentile et al. 2007). This power-law also 
applies, with significant scatter, to elliptical galaxies and bulges if replacing V C i r with ~ 1 — 2 
times the typical stellar dispersion a (Faber & Jackson 1976). Nevertheless, a much tighter 
relation exists for the central black holes of these nearly spherical systems. 

The formation of central black holes (BHs) in galaxies is likely to be a rapid process 
since most quasars have already formed at redshift z > 2. There is a tight correlation 
between the BH mass and the mass of the spheroidal (bulge) component, or even better its 
velocity dispersion (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002). 
The correlations Mb u i ge oc Mbh oc <t 4 are so tight that it is hard to explain unless bulges 
form as fast as BHs, and their growth is controlled simultaneously by some mechanism. At 
the present epoch, the BH accretion rate is both small and completely decoupled from the 
bulge growth. The likely window to couple the two is at high redshift during phases of 
rapid growth and violent feedback. Many previous discussions emphasize that the feedback 
from central supermassive black hole, which interacts with the surrounding environment in a 
self-regulated way, is the key to form the correlations (Silk & Rees 1998; King 2003; Wyithe 
& Loeb 2003; Murray, Quataert & Thompson 2005; Begelman & Nath 2005; Cen 2007). 
On the other hand, the starburst activities peak at similar redshifts to the quasars as a 
whole. In a co-evolution scenario of starbursts and a central BH, the central BH accretes 
with high accretion rate during the main star formation (SF) phase (Alexander 2005). To 
make the starbursts, it was proposed that bulges can form by a rapid collapse due to radial 
instability of isothermal gas. This proposal has the nice feature of forming bulges before 
disks (Xu, Wu, & Zhao 2007). Inspired by these works in a Newtonian framework, we model 
the criteria of bulge formation, assuming a more general mixture of gas and a non-isotropic 
stellar component imbeded in a constant external gravity provided by either CDM halos or 
effective halos of a MOND scalar field. 

Observations show that most of the local and distant starburst galaxies are rich in gas 
and dust (Heckman, Armus, & Miley 1990; Meurer et al. 1995; Sanders & Mirabel 1996; 
Adelberger & Steidel 2000). Photons from newly-formed stars and the BH's accretion disk 
with a luminosity L SF and L BH would diffuse out of the gas sphere. While keeping the 
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gas nearly isothermal, the photons exert a pressure due to dust opacity. The momentum 
deposit rate of photons from an Eddington accreting BH, (2/c)Lbh — (Mbh/ M q ) x 4.3 x 
lO~ 8 ms~ 2 M = (M BH /1.2 x 1O 8 M ) x 10 31 Newton, might drive a feedback force on the 
gas (King & Pounds 2003) to balance the inward momentum deposit of the SF, (2/c)Lsf — 
(Lsf/3.7 x 10 12 L Q ) x 10 31 Newton; the latter could also drive an outward force to balance, 
say half of, the gravitational force on the gas. So we have roughly 

(10' 10 ms- 2 )(430M Bi? ) ~ \ l bh ~ ^L SF ~ 0.5gM g (1) 

where g is the mass-averaged gravity on the gas sphere. Rewriting g = 10~ 10 ms _2 g_io, 
we find a short star formation time scale (0.001c 2 ) M g / L S f ~ 0.2gl{ Gyr if adopting the 
usual SF efficiency of about 0.001 (Leitherer et al. 1999; Bruzual & Chariot 2003). In a 
starburst, we expect radial collapse and violent feedback, so the short timescale of SF should 
be comparable to the free-fall time scale. Assuming that all the gas eventually turned into 
stars, we recover a Magorrian et al. (1998) relation M ^ H ~ 0.002g_io The key point here 
is that observed properties of BH, SF and bulges can all be realised if g-io ~ 2 with little 
scatter universally. 

Presently there are two paradigms where galaxy structure and formation are studied. 
Introducing only two speculative dark components of the universe, the paradigm of Cold Dark 
Matter (CDM) plus the cosmological constant A ~ (10~ 9 m/s 2 ) 2 , it is possible to simulate 
the large scale universe realistically with General Relativity. Many properties of model 
galaxies can be predicted, although not all are in agreement with observations, especially for 
low-surface brightness galaxies. The Modified Newtonian Dynamics (MOND) paradigm is 
able to match properties of high and low surface brightness galaxies with amazing accuracy 
by simply introducing a fundamental scale a ~ 1.2 x 10~ 10 m/s 2 in the space-time metric 
gradient. However, it is generally underdeveloped in terms of simulating the processes of 
structure and galaxy formation. Nevertheless, it was realised that disks in MOND become 
unstable once above certain central surface brightness ~ 13OM pc -2 , so bright regions 
of galaxies tend to exist in spheroidal form and are in strong gravity g > (Sanders & 
McGaugh 2002); throughout the paper we use ao as the dividing scale for strong vs. weak 
gravity. 

The classical MOND gravitational theory (Bekenstein & Milgrom 1984) is now furbished 
with several co-variant versions. These have different constructions using a vector field plus 
(optional) scalar fields (Bekenstein 2004, Sanders 2005, Halle, Zhao, & Li 2008, Zhao & Li 
2008); the scalar field(s) can always be absorbed (e.g., into the modulus of the vector field), 
but can be useful to make the expressions of the theories simpler. In these theories there 
is a constant scale a ~ 1.2 x 10~ 10 m sec -2 ; in regions or epochs of more gradual variation 
of the space-time metric, the dominating source to Einstein's tensor switches from normal 
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matter/radiation fields to the new vector or scalar fields. With the same scale the most 
recent co-variant model (Zhao 2007) can explain how the expansion of the universe switches 
from de-accelerating to accelerating, and explains the amplitude of the cosmological con- 
stant A ~ (8ao) 2 , rather than invoking it arbitrarily as in ACDM. In principle, cosmological 
structure formation all the way to galaxies is well-defined in some of these co-variant the- 
ories (Halle, Zhao & Li 2008). In the original proposal for TeVeS (Bekenstein 2004), there 
is a dis-continuity in the Lagrangian for the scalar field, which makes the evolution from 
cosmological linear perturbations to quasi-static galaxies problematic. This is overcome in 
the proposed Lagrangian of Zhao & Famaey (2006). It is theoretically possible to seed and 
simulate galaxy formation from the linear perturbation simulations of the Cosmic Microwave 
Background (Skordis et al. 2005). Already work has started on galaxy formation modeling 
in MOND (Sanders 2008, Sanders & Land 2008) and there are many constraints on the 
theory from gravitational lensing tests (Feix et al. 2007, 2008, Shan, Feix, Famaey, & Zhao 
2008, Natarajan & Zhao 2008, Tian, Hoekstra, Zhao 2008). 

Here we speculate upon the properties of galaxies formed in the TeVeS framework, 
and contrast with the ACDM framework. While generally speaking the two paradigms are 
mutually exclusive, here we illustrate some similarities, in the context of the formation of high 
surface brightness bulges and elliptical galaxies. There are, however, important differences. 



Let us consider the Cold Dark Matter (CDM) galaxy formation framework. In CDM 
models baryons fall into the potential well of CDM, cool, and condense into stars. The 
background dark matter distribution is often described by the NFW density distribution for 
dark matter (Navarro, Frenk & White 1997, Navarro et al. 2004 ), which decreases smoothly 
from an r~ l density cusp inside to a r~ 3 tail outside a scale radius r s , defined to be where 
the logarithmic density slope is —2 exactly. The density satisfies rp NFW ss p s r s , insensitive 
to radius inside the cusp region r < r s , where p s is a density scale. The Newtonian gravity 
of the halo, determined by the halo mass enclosed inside radius r, is given by 



2. Scatter of the halo gravity in CDM 
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The approximation is valid to 10 percent for < y < 20 and gives F(0) = 1. For a halo of 
virial mass M vir , the halo scale parameters satisfy 

r s = (17kpc/h)M° ir 46 12 , (4) 
r sPs = 13OM pc- 2 xS, (5) 

/ \ 0.34 

where h = H /100 is the dimensionless Hubble constant, and the halo virial mass M vir is 
related to the baryonic mass Mb of a galaxy by M v i r ~ 8M&. The parameters are insensitive 
to the resolution of the simulation and whether the cusp is truly finite or infinite at the 
origin. The numerical values above and the data points shown in Fig. [1] are taken from 
Navarro et al. (1996) and Table 3 of the simulations of Navarro et al. (2004). They agree 
with the scalings for the concentration vs mass, c = ~ 13AM^£\%(1 + z)" 1 , as given in 
Bullock et al. (2001), where c is a ratio betwen the virial radius r„j r and the scale radius r s . 

It is interesting that in the cusp, F ~ 1, and the NFW halo's self-gravity gNFw{r) is 
insensitive to radius and halo parameters, i.e., 

9NFw{r) ~ (1.2 ± 0.6)S X 10- 10 m/s 2 , (7) 

hence gNFw{ r ) of NFW cusps is nearly a universal constant for galaxies of the same baryonic 
total mass. A factor of two scatter in gNFw( r ) is estimated from Table 3 and Fig. la of 
Navarro et al. (2004); the galaxy clusters of 10 14 M Q have nearly the same r~ l extrapolated 
inner density as 1O 12 M galaxy halos (although the cluster simulations do not yet have 
enough resolution to be confidently extrapolated to small radii of a few kpc). 

More precisely r° s m p s ~ est insensitive to halo mass, and simulation resolution, i.e., it 
applies to the r -1,5 Moore profile as well the cored profile of Navarro et al. (2004) as long 
as r s has the model-insensitive definition of the radius where the logarithmic density slope 
is —2. The cluster density is a factor of two higher, and the dwarfs are a factor of four lower 
due to the 5 factor. 

The above is roughly in agreement with Xu et al. (2007), who noted that in the central 
region containing the galaxy bulge, NFW halos have a density scaling relation rp ~ r s p s ~ 
130M Q pc _2 H, and E(M vir , z,c) ~ M^ii ~ 1- While the details differ, in both cases 5 is a 
shallow function of the halo virial mass M V i r , the redshift z, and the concentration c. 
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2.1. Effects on CDM cusp if bulges and ellipticals grow adiabatically 

The CDM density is highly compressible by gravitational interaction with baryons. The 
stellar distribution in elliptical galaxies is often described by Sersic profile in projected light 
(see Graham & Driver 2005), with an underlying volume density of the approximate form 
r -i+o.6/n eX p(_ r 1 /") (Prugniel & Simien 1997); for the Sersic index n ~ 4, we have a nearly 
1/r cusp, which suggests a nearly constant Newtonian gravity in the center. Here we use the 
simpler Hernquist profile for the enclosed mass Mb(r), hence we find the central Newtonian 
gravity has a spatially constant value, 

. . GM b (r) GM b GM b 
9N{r) = 



r 2 (r + Th) 2 r\ ' 

for r — > 0, where M b and are the baryon (total) mass and scale length. The scale length 
and stellar mass (luminosity) of elliptical galaxies are typically correlated with some scatter, 
e.g., Chen & Zhao (2008) find that 

log GU *? - 2 =-1.52 log ( ^±0.1 (9) 

where 1.5 x 10 11 M Q is the characteristic turn-over mass scale in the Schechter stellar mass 

function of galaxies, found by fitting SDSS galaxy samples in the range 10 8 — 1O 12 M 

(Panters et al. 2004), assuming a Hubble parameter 70km/s/Mpc. This would imply 

that near the Hernquist cusp the Newtonian gravity in units of 10 _10 m/s 2 is gN-io = 
( m \~ 1 - 52 

350 x ( i 5 x io" m q ) • ^his mus trate that cores of ellipticals have g ^> a , i.e., they 
strong gravity enviornment (Milgrom & Sanders 2003). @ 



arc 



If in elliptical galaxies baryons condense into the center adiabatically, this process would 
further increase the value of rp^py/ or H. In the most extreme case, a galaxy might start 
with a fb : (1 — /&) = 1 : 7 mix of gas and CDM particles all distributed on circular orbits 
with an NFW distribution of the radius r^, and the gas part condense adiabatically into the 
present stellar mass distribution M b {r). Following the standard recipe (e.g., Klypin, Zhao, 
Somerville 2002), the initial CDM halo mass (1 — fb)MNFw( r i) is squeezed into a radius of 
r, conserving mass. From conservation of angular momentum J of circular orbits, we have 

J 2 = ri GM NFW (ri) = r [(1 - f b )GM NFW (n) + GM b {r)\ . (10) 



1 As an alternative model-insensitive check, we estimate from Faber et al. (1997, their Eq.3-4, and Fig. 4) 
that the typical observed Newtonian gravity near the core or the half-mass radii of cored giant ellipticals 
(with Lio = L/10 10 Lq > 1) is gjv-io ~ IOOL^q ' 6 , and for cusped dwarf ellipticals (with L w < 1) is 
. 9A r_ 10 ~ 100CLL+ ' 6 . 



- 7- 



By Taylor expanding near r ~ 0, we find the contraction factor C = — satisfies 



CT* « (1 - / 6 )C7' + ^- 350 " »1, (11) 

where we considered only the dominant 2nd term and applied the approximation of Chen & 
Zhao (2006). By contracting the same halo mass inside a factor of C smaller radii, one keeps 
the r _1 profile of of the halo cusp, but increases the halo central density and self-gravity by 
a factor C 2 , i.e., at the center 

9cdm = 2-kGtpcdm = C 2 g NFW = C 2 E ■ lQ- w m/s 2 > a . (12) 

The above is likely an over-estimate, since more realistic formation of ellipticals will involve 
mergers of gas-rich spirals galaxies. In any case, there is likely a large upward scatter around 
the naive analytical prediction gcDM ~ 10~ 10 m/s 2 due to different scenarios of formation 
and different halo masses and concentrations etc.. 



3. A universal gravity scale: effective halos in simple MOND 

In the co-variant theories of MOND, the only sources of gravity are the stars and the 
gas (i.e., the baryons). In the context of its co- variant version, the Tensor- Vector-Scalar 
theory of Bekenstein (2004), there should be a scalar field <f> a , such that in the spherical case 
g s = — V0 S gives a dark halo like gravity. Here we call the scalar field the Effective Dark 
Matter (EDM), qedm = 9 s = |gs|; it is related to the Newtonian gravity g n = |gjv| an d the 
actual acceleration (or gravity) g = |g| = + g s |. The Poisson equation is modified as 

V ■ (tJL 8& ,) = V ■ (fig) = V ■ g N = -4nGp, (13) 

where /i s and fi are modification functions, which give identical descriptions in case of a 
spherical mass distribution, and reduce to Newtonian dynamics when g^ S> ao- 

Consider the modification functions as proposed in Angus et al. (2006), 

9s _ CLq 

li s = — , a a = — , (14) 

{a a - 9s)ct a 

where a a is the fixed scale of the theory with a being a theory constant; the a = 1 "simple" 
model is the most popular special case. 

Reexpressing g s in terms of a rescaled Newtonian baryonic gravity y, we find that § 



2 Here /i = 1 



g+ag 



n -1 



/ g-a a \ I g 

\ 2a a J a.a a 



The combined gravity of effective DM and baryonic matter 



is the actual acceleration ^ = g(r) = a a [0 s (y) + a 1 y\, where s (y) = ^— j=g _ 1 ■ Note 9 s (y) ~ 1 if y 3> 1. 
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around a gas plus stellar sphere there will be an effective DM halo gravity Qedm or the 
scalar field g s , 

G(M g + M> 

9edm = 9s = a a 9 s {y), V = 2 (15) 

where M g + M* is the gas plus star mass inside radius r. A remarkable result of the general 
class of /i-function is that there is a maximum to the scalar field gravity 

9s < 9s,max = a a . (16) 

This means that when the Newtonian gravity, g^, is the strongest (as in the centers of bright 
galaxies), the scalar field g s — > a /a, i.e., approaching a universal constant plateau. This 
breaks down only if a = 0, corresponding to Bekenstein's toy function. Also if a — > oo, 
the dynamics becomes purely Newtonian with a zero scalar field. The standard fi(x) = 
x j \J\ + x 2 1 which fits observations well, can be approximated by a ~ 3 in terms of sharpness 
of transition from strong to weak gravity. One can set a ~ 1 — 3 to be consistent with galaxy 
data. | 

The density profile of the EDM can be derived from the Poisson equation as 

PedmM = ^(r 2 g EDM ), . (17) 

Taking gEDM = a a @s and the Newtonian gravity of the baryonic Hernquist profile (Eq. EJ), 
we find 

_ da „ a -fa d9 ° r 
PEDM - 7, TT^U Vl = [Va - -J] — 

ZiraGr \ amy r + T} lj 

For small radii (~ kpc) Pedm has a roughly 1/r cusp because Q\ ~ 6 S ~ 1 at centers of 
bright galaxies, where g N ~ GMjr\ ~ 350(M/1.5 x lO n M )~ L52 a o > a . In fact near the 
center g N 3> a ~ a a even if adopting a scale length 3 times bigger than implied by the mean 
scaling in Chen & Zhao (2006), and/or using a Sersic profile for the stellar distribution. The 
strong gravity implies a saturated scalar field in the center: g s reaches its maximum a a . This 
is confirmed numerically at least for the a = 1 case (cf. Fig. [2]). In fact near the centers of 
ellipticals the model predicts 

T Pedm = II = ~ 13OM p C - 2 . (19) 



3 Based on theoretical arguments and matching observed rotation curves, Zhao & Famaey (2006) advo- 
cated the "simple" function with a = 1, fi s = ao ^ g , [i = ao 9 +g ■ This /^-function is also supported by 
the Milky Way kinematics data, and by the SDSS extragalactic satellite velocity distribution (Angus et al. 
2007), and is consistent with the recent rotation curve fittings by Famaey et al. (2007a, b), Zhao & Famaey 
(2006), and Sanders & Noordermeer (2007), McGaugh (2008). 
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This is rather similiar to the case of NFW halos, but rpEDM has virtually no scatter for the 
MONDian scalar field in bright centers of ellipticals. Fig. [2] shows the scalar field is very 
rig 

id, B g s = (0.5 - l)a /a, nearly incompressible on scales of 1-10 kpc, i.e., it cannot be 
increased significantly by compressing the baryonic material and increasing the Newtonian 
gravity. We can also define a maximum central pressure of the scalar field for later use [f| 

a 2 

P Q = a a = ao/a. (20) 

It is remarkable that the centers of ellipticals are all immersed in strong gravity, hence 
a universal constant scalar field, perhaps extending all the way to the central black holes of 
these system. One ponders the consequences of such a remarkable uniformity and univer- 
sality for central kpc regions of all bright galaxies. One wonders, in particular, whether the 
uniformity would lead to very tight correlations, such as the Mbh — o~ relation. 



4. Modeling spheroid formation in a constant background gravity a a 

As a demonstration of the consequence of a universal scale, we apply it to galaxy for- 
mation and scaling relations. Unless stated otherwise we make the approximation 

g NF w ~ 9edm = 9s = a a 6 s ~ a a = est (21) 

This approximation means that baryons experience a rigid uniform extra field a a on top of 
the Newtonian self-gravity of the stars and gas, 

g(r)=a a + ^^- 2 *->-. (22) 

• If the origin of the extra field is actually from NFW cusped CDM, then it should be 
understood that a a would have a significant scatter between galaxies with a trend for 
bigger a a for bigger galaxies. 

• If the origin is the scalar field in a co-variant version of MOND, then a a should be 
viewed as a universal fundamental constant intrinsic to a theory. It is equal to a Q /a 



4 The scalar field becomes half-saturated as soon as the overall gravity exceeds ao/a, where g s — 9n = 
0.5ao/a. In the Milky Way this translates to a galactocentric radius of about (220kms -1 ) 2 a/(ao) ~ 13akpc, 
which is the edge of the galaxy disk. 

5 Zhao (2007) noted that this pressure P a ultimately relates to the cosmological constant or the vacuum 
energy density, which has the same order of magnitude as P a 
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with zero scatter, although the value of a ~ (1 — 3) is not precisely determined at 
present. 

Spatial non-uniformness of the background field is studied in the Appendix; the effect can 
be treated crudely as a small spatial variation of a. 



4.1. Maximum stable gas mass 



First we construct analytical spherical models of gas and star equilibrium in MOND, 
and study the condition for the gas sphere to remain stable. It is found that as we increase 
the core density up to a critical value, the total gas mass increases. Once the gas core density 
exceeds the critical, the total gas mass starts to decrease. 

Assume a quasi-static equilibrium of a gas sphere p g (r) of sound speed c g {r) and a stellar 

2 

sphere p*(r) of radial velocity dispersion cr* and an anisotropy parameter (3=1 ^. To 

model the tangential velocity dispersion and non-isothermal profile we introduce a velocity 
dispersion measure <7i, 



2 _d\ogp*al d\ogr 

0-1 = 4*0-*, 4* = — — (J. 

a log p* a log p* 



(23) 



Define £ to be the ratio of thermal pressure p g o~g to the opacity-induced pressure (p s cr|) / £ 
on the dusty gas sphere (not the stellar sphere), countering the gravity. The velocity dis- 



persion o-gij-) is related to the sound speed by c 
mixture satisfies the equations 



9 d log p g 



In equilibrium the gas-stars 



9(r) 



_£ 

r 



-d\n[p g (r)\ 
d\nr 
dln[p*(r)] 



dM g (r) 



d\nr 

dMJr 



;i+r 1 ) 



(24) 



p g (r)dr p*(r)dr' 



(25) 



The conversion of gas into stars also needs to be modeled to be realistic. The simplest 
solution of the above eqs. would be a model where stars trace the gas radial distribution, so 
we have 

p*(r) M*(r) /*(*) 



M s (r) 



!-/*(*) 



(26) 
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where the position- independent factor /*(£) is the fraction of gas formed into stars at time 
t. Such a solution is possible if the feedback ratio £ _1 is regulated by star formation, and 
related to the temperature of the gas by 

(1 + r 1 )^ = = o\ = est. (27) 

Note that the feedback parameter £ and the gas sound speed c g are not required to be 
rigorously independent of radius, although certain combination of these two is. We do not 
require the stellar component to be exactly isothermal or isotropic either. 

We compute the gas equilibrium for a range of core pressures p(0) (see Figure [3]) after 
rewriting the equations in term of the dimensionless mass m, radius u(m) and rescaled 
density p(m) (see Appendix). We find that the gas density generally falls off monotonically 
with radius or mass. All models have finite mass out to infinite radius. The density drops 
steeply with radius due to the deep linear potential well of the background gravity, hence the 
mass converges quickly. Curiously there is also a maximum m max « 4.3 in the total rescaled 
mass as p(0) increases. This happens at a critical core density or pressure 

p( ) = Ml « 30 (28) 

above which the gas density p g of a parcel of gas dM g no longer increases monotonically with 
an increase of the central pressure, and in fact the total mass will decrease with increasing 
p(0) after it reaches the maximum value. 

These limits on gas central pressure and total mass are examples of the instability first 
discussed by Elmegreen (1999). A gas sphere above a certain critical mass 

M max « 4.3 f-p^j , (29) 

or above a critical central gas density or pressure, does not have stable solutions; adding a 
tiny amount of gas would lead to collapse. It is interesting to speculate that bulge formation 
originates from such a gas instability. 

For a sphere of gas plus stars at the critical mass, we integrate the density numerically, 
and find the central surface mass density 

5(0) ~ 2a a /G ~ lSOOa^Mopc^ 2 , (30) 

insensitive to the initial gas velocity dispersion a g and the feedback ratio We also fit 
our numerical solution of the gas mass profile by a Sersic-like distribution. The gas mass 
profile is related to the dimensionless volume density profile j(u) via 

MJr,t) [ u , . . „ , , r 

f J [ = / j(u(Wdri, ^ =-!-!■ (31) 

(1 - }*{t))M max J a(a a L 
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Prugniel & Simien (1997, eq B6) suggested an approximated (valid to 5% accuracy typically) 
form r - 1 +°- 6 /™ eX p(— r 1 /™) for the deprojected volume density of a Sersic profile. We fit j(u) 
numerically to such a Sersic profile of volume density, and find an — 1.2 profile works well, 
i.e., 

j(u) « 0.14w~ 1/2 exp(-lW /L2 ). (32) 
The normalisation is such that J °° 4irj(u)u 2 du = 1. 



4.2. Properties of the stellar component formed 

In the simplest scenario gas might turn adiabatically into stars while maintaining the 
density profile at the critical density and mass. Eventually /* = 1 when the gas is exhausted 
by star formation, we have a stellar system with a nearly Sersic profile of high central surface 
brightness. The density slope is shallower /steepr than isothermal inside/outside the radius 
r ~ O.Safa- 1 . 

Our model also resembles real spheroids since they satisfy the Faber- Jackson-like relation 
M£° ~ between the total mass and stellar velocity dispersion (Faber & Jackson 1976), 

M " = ^ = 4 - 3 ^ xl ° n ^(») 4 - (33 > 

Considering small deviations from isothermal and isotropic f3 ~ stellar distribution, cr* = 
0"i/£* can be treated as effectively certain mean projected stellar dispersion. Note that our 
result differs in detail from the MOND virial theorem M* = 4^-0"* (Sanders & McGaugh 
2002) which applies to an isotropic isothermal stellar system in deep-MOND (see Appendix). 
Our bulges and ellipticals are clearly in a mild or high acceleration regime. 



4.3. Correlations of Black Hole, Stellar Spheroid and Starburst 

Similar to the introduction, we assume the momentum deposit rate of photons, -(Lsf + 
Lbh) drives a force acting on the gas. Assume the BH grows exponentially by Eddington 
accretion from a seed before reaching the balance Lbh — L$f- After this point, the BH 
mass stops growing and we set Lbh = 0, and the gas mass and the Lsf will then exhaust 
exponentially as the starburst continues. Assuming at the turning point, the BH and SF 
each provide half of the total feedback force, which is a fraction (1 + ~ 1 of the radially- 
pointing gravity g(r), we obtain 

\ L ™ = l LsF = Wh) m ' (34) 
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where M g = (1 — f*)M max is the total gas mass, and the mass-averaged gravity 

g = M; 1 J gdM g ~2a a (35) 

is computed by an integration of mass for the density model with the critical mass. This 
result is robust for a range of the polytropic index of the gas (see Appendix). 

Observations reveal that some gas rich local spirals have /* ~ 0.5 ± 0.2. The stellar 
fraction is /* ~ 1 for local bulges and elliptical galaxies. This is much higher than in their 
z ~ 3 progenitor starburst galaxies /* ~ 0.02 (Alexander et al. 2005; Genzel et al. 2006; 
Kennicutt 1998, Solomon & Vanden Bout 2005). We shall adopt M g = (1 - f*)M max = 
(0.5 - 1)M~, and g ~ 2a a ~ (1 - 2) x lO^ms" 2 , and (1 + £)~ X ~ 1 as long as the feedback 
l/£ is not too small (Murray, Quartet & Thompson 2005). We obtain in the end 

Mbh L sf _ M?? ( C*. \ 4 

2 x 1O 8 M 6 x 1O 12 L a x 10 n M o V 200 kms_1 / ' 

where £' — (1 — /*)/(l + £) ~ 0.5 — 1, consistent with the observed Magorrian relation 
M BH ~ 0.005M bulgc , and = (4.3£') 1/4 £* ~ 1. The scatter in ^ 4 is about factor of a 
few each way, consistent with the narrow scatter seen in the observed relation of Mbh ~ 
(1 -4) x lO 8 M (o-,/2OOkms~ 1 ) 4 (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine 
et al. 2002). 



4.4. Hydro simulations of gas collapse in a fixed background potential 

To confirm the analytical results on the gas instability, we investigate the threshold of 
gas collapse numerically. We consider the gravitational collapse of gas within a rigid dark 
matter halo using a hydrodynamic simulation. The gas which forms the bulge is embedded 
in a uniform external field from the dark matter potential. The simulation represents a 
simple treatment of the CDM halo gravity or the MOND scalar field. 

This simulation was performed using the Lagrangian fluids code, Smoothed Particle 
Hydrodynamics (SPH). The code is described in Bate, Bonnell & Price (1995), and is based 
on the version by Benz et. al. (1990). The simulation uses 10 5 particles, which are initially 
distributed randomly within a sphere of radius 15 kpc. The gas is initially at rest, with a 
temperature of 10 6 K. The mean molecular weight is 1.2, thus the sound speed of the gas 
is 83 km s -1 . The dark matter halo is included by way of a fixed linear external potential. 
The potential per unit mass is — a a r such that there is a universal acceleration of a a = 10 10 
m sec -2 towards the centre of the galaxy. The gas is subject to both the external potential 
and self gravity during the calculations. 
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With these parameters, the critical mass is 1.54 x 10 10 M . We first set the total mass 
of the sphere to 0.0065 M max (10 8 M ), so the gas should be stable against gravitational 
collapse. We then evolved the gas for 12 crossing times, during which time the gas settles 
into equilibrium according to the dark matter potential. The mass was then doubled (by 
doubling the mass of each particle) and the calculation resumed. By 5 crossing times, the 
gas had again settled into equilibrium. This process was repeated, doubling the mass each 
time the gas reaches equilibrium, until runaway gravitational collapse occurs. Fig. H] shows 
snapshots of the radial profile of the gas pressure just below (upper panles) and just above 
(lower panels) the critical mass. A collapse of gas is evident from the density cusp in the 
final state. 

In Fig [3J a 1 D projection of the column density (along the x axis) is plotted versus 
radius. During the first run, the distribution of gas changes from a uniform profile, to a profile 
corresponding to the dark matter halo. Over the subsequent runs, there is little change in 
the distribution whilst the total mass remains less than the critical mass. However when the 
mass exceeds the critical mass, the cloud is gravitationally unstable. Gravitational collapse 
occurs and the density at the centre of the cloud continuously increases. The calculation 
stops, the final profile shown as the thickest line in Fig. [5J the slightly subcritical gas can be 
approximated by a Sersic law with n = 1 (i.e., exponential). 

Galaxy formation clearly involves more than just spherical collapse. The rapid spherical 
collapse phase might produce the bulge and the black hole. This phase is likely followed by 
a gradual phase of episodes of minor mergers, the dense galaxy will likely acquire a diffuse 
stellar halo of higher angular momentum material. One might expect a shallower outer 
profile than n = 1 in the end depending on the amount of stars accreted. Indeed real 
galaxies have a range in Sersic index, which is correlated with the stellar mass of the galaxy 
(Caon et al. 1993, Desroches et al. 2007). Bulges and pesudo-bulges typically have profiles 
with 4 > n > 1, and the most massive elliptical galaxies have profiles with n > 4. Some 
discrepancy with observations is expected since a constant scalar field is a poor assumption 
at large radii. 

5. Summary 

In short, combining the results of analytical models and numerical simulations, we argue 
that gas collapse and BH feedback inside the effective DM halos of MOND can produce 
galaxies with realistic scaling relations. Better than in Newtonian NFW halos, the Faber- 
Jackson and BH mass-velocity dispersion scaling relations are recovered with narrow scatter 
thanks to the fact that there is a redshift-insensitive and luminosity-insensitive universal 
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scale of gravity in the high-z gas-rich starburst galaxies and present day ellipticals. Many 
feedback processes (e.g., Silk & Rees 1998, King & Pounds 2003, Cen 2007, Xu, Wu, Zhao 
2007) invoked in the CDM context could be important in MOND context as well, and could 
change predictions of the state of the art MOND N-body simulation (Nipoti et al. 2008) 
and hydro simulations (Tiret & Combes 2008). 

In the context of a simplified picture of TeVeS, a nearly-isothermal gas sphere can 
collapse and trigger a starburst if the gas central pressure is above a universal threshold. 
This condition is likely synchronised throughout the universe, consistent with the observed 
epoch of starbursts. We also recover the Mbh — o"* relation, if the gas collapse is regulated 
or resisted by the feedback from radiation from the central BH. 

While the pristine CDM halos give acceleration of typically the order (1 — 3)a , there 
are two important differences with the universal scale in MOND. In CDM, this scale would 
exhibit a large variation due to scatter in halo concentration (Milgrom 2002). Secondly, as 
the galaxy bulge forms, the central part becomes baryon-dominated in Newtonian gravity 
and the CDM is adiabatically compressed to higher densities without a theoretical upper 
limit unless DM is made of neutrino or its sterile partners (cf. Angus 2008, Zhao 2008). 
The collapse threshold could rise as the background gravity a a increases. In contrast, in 
the TeVeS picture the scalar field g s within all bright galaxies stays close to a a throughout 
galaxy formation. This value universal constant in MOND. Hence we expect a 

universal theshold and synchronised formation of galaxy spheroids at high redshift, perhaps 
producing star burst galaxies in TeVeS. The mechanism would work less well in the CDM 
paradigm. 

Speculating beyond our models, elliptical galaxies are often triaxial in their stellar dis- 
tribution and potential. If collapse happened above the threshold described here, the uneven 
distribution of angular momentum would lead to a triaxial equilibrium potential. This expec- 
tation is consistent with N-body simulations of collapse in CDM halos (Navarro et al. 1996) 
and in MOND (Nipoti et al. 2007). Indeed, elliptical galaxies with a mild r" 1 cusp in stellar 
light are allowed to exist in self-consistent triaxial configuration in CDM (Capuzzo-Dolcetta 
et al. 2007) and in MOND (Wang et al. 2008). 

Corrections are expected for low surface brightness galaxies since the scalar field in these 
systems are far from being saturated. So g ~ wa a , where w < 0.3, and we expect the faint 
dwarf spheroidal mass and their central BH mass to reduce by a factor w and w 2 respectively 
compared to a blind application of Faber- Jackson and M BH — a relations. No correction is 
expected for M32-like compact dwarf ellipticals, which are observed to have bright stellar 
nuclei and BHs. Our prediction that the BH mass per unit stellar luminosity Mbh/L* is 
lower in dwarf spheroidals than dwarf ellipticals could be checked by sensitive searches for 
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(central tracers of) massive BHs in dwarf spheroidals, which is challenging observationally. 

Environmental effects could lead to corrections to the Faber- Jackson relation. Member 
galaxies in the center of a rich galaxy cluster have orbital accelerations of ~ 2ao, which makes 
it easier to saturate their MOND scalar field due to the external field effect; any gas rich 
system would be more vulnerable to radial collapse into a high surface brightness triaxial 
galaxy in the center of a galaxy cluster than in the surrounding or in the field (Wu et al. 
2007, 2008). The external field effect is similar to making a a smaller. Interestingly, eq. (1301) 
and f)33p can be combined to form 

Ml <_^A ( a * \Y 1500L qp c " 2 V- 3 ^ 1q11l q rs7i 

~ hL* G 2 ~ \200kms~ 1 J \ /* J L* 1 ' 

where i* = S*/(M/L) is the surface brightness of stellar distribution of surface density S*. 
This relation is very similar to the fundamental plane relation of ellipticals (e.g. Binney & 
Merrifield 1998) 

4\ 0.7 



L* oc [f\ , (38) 

if we accept that M*/L* ~ L° 2 . Equivalently we can write the galaxy size r e oc S^V 2 , which 
is very close to the lensing mass fundamental plane oc <J 2 /r e determined by GR-based 
strong lensing results, S® m oc a im /r e (Bolton et al. 2007); the MONDian corrections to 
the lensing mass within the Einstein radius are generally mild (Zhao, Bacon, Taylor, Home 
2006, Shan, Feix, Famaey, Zhao 2008). On the other hand, environmental effects would not 
correct the BH-velocity dispersion relation because 

M BH oc at, (39) 

which is independent of the parameter a a . 

The BH-bulge relation may have been frozen at high redshift; some scatter might be 
caused by gas rich mergers, and minor feeding of the BH through a rotating gas-rich bar later 
on, which might evolve secularly into a bulge (Tiret & Combes 2007). Disk galaxies could 
have built up their disks from high-angular momentum gas well after the bright bulges form 
by radial collapse. The formation of the disk part will not fuel the central BH significantly. 
This may explain why the BH mass is tightly related to the bulge and its velocity dispersion 
rather than with the total baryonic mass and the terminal circular velocity in disk galaxies; 
the latter two are correlated themselves by the Tully-Fisher relation, which is built in the 
Lagrangian of covariant MOND because virtually no scatter from this relation is observed 
in field galaxies (cf. Wu et al. 2007). 
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Fig. 1. — Shows the scatter in logr s in kpc vs log p s r s in M Q pc -2 in simulated NFW halos. 
The line shows our empirical relation p s r s = 130HM Q pc~ 2 , and S = 2(r s /17/i~ 1 kpc) a34 (cf. 
eq.5). The data points from Navarro et al. (2004, adopting = 0.7 and a 8 = 0.9, h = 0.7) 
are systematically below those from the original NFW (1996, adopting = 0, as = 0.63 
and h = 0.5); apart from scatter, the halo concentration has systematic dependence both 
on the cosmology and on the halo size. In comparison MOND predicts a universal effective 
DM scale log(rp) ~ log 130 ~ 2.3 (not shown). 
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Fig. 2. — The mass distribution of Effective DM halos (left panel and assuming an a = 
1 MOND) and NFW Cold DM halos (right panel) in terms of rp(r). Note the nearly 
universal central density, especially in MOND. Not shown here is a further factor of 10 scatter 
intrinsic to NFW CDM halo density, and another factor of ~ 10 — 100 upward correction if 
elliptical galaxies form adiabatically and compress CDM. For MONDian rp EDM) four lines 
represent M b = (10 12 , 10 11 , 10 10 )M Q , respectively, as indicated. For CDM rp NFW , we take 
the corresponding value of the halo mass as M NF w = 8M b . 
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Fig. 3. — Mass distribution of a gas sphere embedded in a rigid background uniform gravity 
a a for models (from bottom to top at small radii) with increasing dimensionless central gas 
pressure p = 1, 3, 10 (bottom three), and p = 30, 100 (top two). The axes are logp(m) vs 
the logarithm of the rescaled radius u(m) = a y a , where m = Ga a o^M g is the rescaled 

baryonic mass, pirn) = ff/^G) is the rescaled gas density or pressure. Note how the high-p 
lines are above/below the lower po lines at small/large radii respectively. This reversal is a 
feature of reaching a maximum in total gas mass at the critical pressure (p ~ 30). Also 
shown is a Sersic {n = 1.2) profile (diamonds). 




5 10 5 10 

radius radius 

Fig. 4. — shows the gas pressure radial profile of the 3D hydro simulations in a fixed back- 
ground field a a = a . Panels (a) and (b) are the "before" and "after" evolution of a gas 
sphere of low mass M g = 3.57 in units of af/(Ga a ) (see text for physical units), which is 
below the critical value 4.3. Panel (c) and (d) show the prominent development of a gas 
density cusp for a high mass M g = 7.14, exceeding the critical gas mass. Note the gas 
profiles are nearly exponential and are consistent with a Sersic index n about 1-1.2 before 
it becomes unstable. In d), the pressure in the centre exceeds the critical pressure, which is 
~ 3.6 x 1(T 15 bar. 
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Fig. 5. — Shows hydro simulation in a fixed background field a a = a . The numbers on the 
left are total gas mass M gas in units of af/(Ga a ). The dashed lines represent the input gas 
profile 'before' the evolution, and solid 'after' the hydro evolution. The timescale for each 
evolution is 5 crossing times, except for the first run from a uniform distribution, which was 
evolved for 12 crossing times. The radius is in unit of Rq = af/a a , and the projected density 
in units of So = a a /(2irG). Instability happens when the critical mass exceeds 4.3 in these 
units, or the surface density exceeds about 10 in these units. The slightly subcritical gas can 
be approximated by a Sersic surface density with n = 1 (plus signs). 
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Fig. 6. — Effective dark matter (scalar field) gravity g s in units of MOND a a (a = 1) as a 
function of enclosed gas mass in units of a g (0) A /(Ga a ) for a non-isothermal polytropic gas 
sphere of central velocity dispersion cr(0) and with 7 = 1.1 (rightmost curve), 1.2, 1.33, 1.5, 
1.66 (leftmost curve). All the curves have central pressure 30P a , and the pressure drops to 
zero at the last points of the curves (truncation radii of the polytrope). 
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A. Mean gravity and scalar field in MONDian polytropic gas models 

with/ without feedback 

A.l. Constant background gravity model with Feedback 

We rewrite the hydrostatic equations using the dimensionless quantities 

_ M g + M* _ r _ p g ( r )+p,(r) 

and expressing the rescaled mass m as the independent coordinate, the problem is recast to 
solving the pair of dimensionless ODEs 

u 2 dp(m) „ m u 2 du(m) 1 ... 

7^ = e > + "7 — \2> a = -r^- A2 

dm u[m) z dm p[m) 

For the constant background gravity models, we set 6 S = 1. For each value of p(0), the 
density profile under the hydrostatic equilibrium can be completely determined with the 
initial conditions at the center for the radius w(0) = and the rescaled density p(0). Results 
are shown in Fig. [21 



A. 2. Polytropic self-gravitating models 

The above simplified models apply only if we can treat the scalar field as constant, and 
the velocity dispersion as isothermal. These are not rigorous for real gas. To check the 
robustness of our results we have also computed numerically more realistic models where we 
solve for a self-gravitating polytrope ^ 2 g p l g ~^ = est in rigorous self-consistent equilibrium in 
MOND gravity. There are no SF and BH in these models. Nevertheless the pressure force 
on the gas sphere 



o 



F = I ( P y g )d(^ r ) (A3) 

\A^)d{p g a 2 g ) (A4) 

(4irr 2 )p g gdr (A5) 

g(r)dM g (A6) 
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where we applied integration by parts, and the hydrostatic equilibrium equation 

d{p g a 2 



dr 



9P 9 - (A7) 



The radial distribution is modelled in the style of Sanders (2000), but here with the 
simple scalar field interpolation function fi s (g s ) = g s /{ a o ~ otg s ), so 



m „ „ am 
+ S 



u(m) 2 u(m) 2 ' 



(A8) 



where m and u(m) are the dimensionless mass and radius, 



GM(r)a a . ra a 

m = , u[m) = — , (A9) 

a\ erf 

where o\ is a characteristic dispersion of the polytrope. 

For models with a — 1, the total gravity g = (1 + fi s )g s = a °_ gs , we computed (see Fig. 
[6]) the scalar field strength for models with a central pressure ~ 30P Q . These models have a 
finite radial truncation, where the polytrope temperature <j 2 g falls to zero. 

The profiles vary with the polytropic index 7, a value typically between 1 and 5/3 
for real gas. Nevertheless the scalar field is roughly constant even in these more realistic 
models where the velocity dispersion is not isothermal. Numerical integration finds that the 
mass-averaged gravity 

9 = Jff 9(r)dM g ~ 2a Q , (A10) 

where we used the approximation that at half-mass radii, g s ~ 0.66a a , = fi s g s ~ 2g s , 
and [i = j^ 3 — ~ 0.66. So our simplistic treatment should be a good guide to real galaxies 
in MOND; extending our self-consistent hydrodynamics code for Fig. H] to incorporate a live 
MOND scalar field would allow us to check the universality of the threshold here. 



A. 3. Consistency with earlier models 

As a consistency check with earlier MOND models, we note that our finding of a maxi- 
mum scalar field acceleration 

\g - 9n\ < a a = — (All) 

a 



6 Alternative expressions are F = gM g — J ^f^dMg — J 2cTg ^ dM g , which applies to non-isothermal gas 
as well. 
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is consistent with Brada and Milgrom (1994) 's first finding of a maximum effective halo 
acceleration in many MOND /^-functions, especially the standard fi function. Fig. [2] argues 
the scalar field is largely a constant in the bulge regions (1-10 kpc). We noted some similarity 
of a NFW halo to the MONDian scalar field with the /i-function of Angus et al. (2006), 
especially if a = 1. The scalar field in the TeVeS picture can be simply replaced by g s — > 
|y9o|^|_£M j n Q^gj. co- variant metric versions of MOND (e.g., Zhao 2007); the counterpart 
in the modified inertia picture is less clear. Earlier authors have also noted the similarity 
of the standard /i-function with Pseudoisothermal halos (Sanders and Milgrom 2005), where 
the halo acceleration goes up to a maximum as a solid body, and then drops. Several 
literatures claim that NFW halos produce somewhat higher than ao accelerations (e.g. Fig.8 
of McGaugh 2004), but the comparison of NFW and MOND has been as systematic as done 
here. Works of Famaey and Binney (2004) and Sanders and Noordmeer (2007) show that 
real data tolerates the a = 1 /i-function. Using the standard /i, Milgrom (1984) and Sanders 
(2000) also note a maximum mass of an isothermal self-consistent distribution in MOND, 

4 4 

about (15 — 20)^7 ~ 4.3a-p^ if we approximate a ~ 3. Their assumption is different 
because their scalar field is not fixed as a constant, but has a peak value about a /a ~ a /3. 
While differing in details, these earlier analysis are qualitatively consistent with our finding 
of an maximum scalar field of order a^/a in MOND. 



